Scaling in a cellular automaton model of earthquake faults 
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We present theoretical arguments and simulation data indicating that the scaling of earthquake 
events in models of faults with long-range stress transfer is composed of at least three distinct 
regions. These regions correspond to three classes of earthquakes with different underlying physical 
mechanisms. In addition to the events that exhibit scaling, there are larger "breakout" events that 
are not on the scaling plot. We discuss the interpretation of these events as fluctuations in the 
vicinity of a spinodal critical point. 
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Earthquake faults and fault systems are known to ex- 
hibit scaling jl]j|] where the number Nm of earthquakes 
with seismic moment M scales as Nm ~ 1/M B with B 
between 1.5 and 2.0 ||. The observed scaling is over sev- 
eral decades, but for the larger events there is an indica- 
tion that scaling does not apply, a fact often attributed 
to poor statistics. However, because models also pro- 
duce this deviation from scaling, even when there are 
many large events the origin of this deviation lies 

elsewhere. Other questions of interest include: What is 
the physical mechanism that produces the scaling? Do all 
the events on the scaling plots have the same physical ori- 
gin? We report the results of our theoretical and numer- 
ical investigations of a cellular automaton (CA) model of 
an earthquake fault indicating that the scaling region is 
dominated by a spinodal-like (pseudospinodal) singular- 
ity that determines the distribution of events. The scal- 
ing can be decomposed into three distinct regions driven 
by different physical mechanisms. In addition to the scal- 
ing region, we find that the largest "earthquakes" are not 
on the scaling plot and have yet another physical origin. 

The system of interest is a CA version of the slider 
block model j| and consists of a discrete two-dimensional 
(d = 2) array of blocks connected by linear springs with 
a spring constant (stress Green's function) T(ry-) and to 
a loader plate by linear springs with constant Kr,; rij 
is the distance between blocks. Each block i initially 
receives a random position Ui from a uniform distribu- 
tion, and the loader plate contribution to the stress is 
set to 0. The stress <Tj on each block is given by a{(t) = 
Ej T{no)Pi (*) ~ W)\ + K l [V £„ 6(n- 1) - U t (t)] and 
compared to a threshold value af. If Oi < erf , the block 
is not moved. If Oi > af, the block slips (fails) and is 
moved according to Ui{t + 1) = U l {t) + [a t {t) - af{t)]/K, 
where K = K L +K C , and K c = E 3 ^ 3 T ( r v)- The 
residual stress, af(t) = a R + a(r]i(t) — 0.5), specifies the 
stress on a block immediately after failure. The random 
noise r\i is taken from a uniform distribution between 
and 1, a sets the noise amplitude, and a R is the average 
residual stress. After all the blocks have been tested and 



moved, the stress on each block is measured again and 
the process is repeated. We choose — Kc/q for all 
j inside a square interaction range with area (2R + l) 2 
centered on site i, where q = (2R + l) 2 — 1 is the number 
of neighbors; = for all the sites outside the interac- 
tion range. After block i slips, Kq/K of the local stress 
drop, <Ji — erf', is distributed equally to its neighbors, and 
Kl/K is dissipated. When no block has a stress greater 
than af , the earthquake ceases and the seismic moment 
released during the event is M = ^ A[/„ where At/, 
is the slip of block i during the earthquake. The loader 
plate is then moved a distance VAT, the stresses are up- 
dated, and we search for the unstable blocks that will 
initiate the next event. The quantity AT, which we set 
equal to 1, sets the "tectonic" time scale. In the limit 
V = the stress is globally incremented to bring the 
"weakest" block to failure and there is a single initiator 
per plate update. 

Because the Ty appropriate for earthquake faults is 
long-range ||, we will consider R » 1. In our sim- 
ulations R = 30, af = o F = 1 is a spatial constant, 
Kl = 1, Kc = 100, V = 0, and the distribution of resid- 
ual stresses is defined by a R = 0.25 and a = 0.5. In Fig. | 
we plot the log (base 10) of the probability n(s) of events 
of size s (number of failing blocks) versus log(s) gener- 
ated by the model. For the chosen parameters there are 
no multiple failures of the same block during an earth- 
quake and M ~ s. For the total of 18 x 10 6 events, there 
is still a significant spread of the data in the large events 
region. The origin of this spread is not poor statistics. 

We now review the theoretical arguments that describe 
the scaling of events. In the limit that R diverges such 
that Jr 2 T(r)dr — > oo but fT(r)dr is finite, we have 
derived a Langevin equation for the stress This 
derivation and numerical simulations |^^] confirm that 
the CA model is described by equilibrium statistical me- 
chanics in the limit of R — > oo and that this description is 
a very good approximation for systems with long, but not 
infinite, R. Because this Langevin equation is a general 
description of systems with a simple scalar order param- 
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FIG. 1. Log-log plot of the probability n(s), the number 
of events with s failing blocks divided by their total number, 
18 x 10 6 . The system consists of 256 x 256 blocks with periodic 
boundary conditions. Note the deviation from the straight 
line of slope 1.5 ± 0.05 for large s. 



eter j8|, the scaling of the fluctuations in the vicinity of 
spinodals of mean-field Ising models (and simple fluids) 
and the present CA model is the same. 

Our main assumption is that the structure and dynam- 
ics of earthquake events is identical to the structure and 
dynamics of fluctuations near critical points and spin- 
odals. In particular, scaling is determined by the pres- 
ence of a spinodal singularity (|||. For Ising systems, by 
mapping the thermal critical point onto a properly cho- 
sen percolation model |||l0| , the properties of the fluctu- 
ations at the thermal critical point can be obtained from 
the properties of the clusters at the percolation thresh- 
old: percolation clusters are the physical realization of 
the fluctuations PJlO[|. At the critical point the clus- 
ters associated with the divergent connectedness length 
are the fluctuations associated with the divergent sus- 
ceptibility in the thermal model. We can use a similar 
mapping to generate a percolation model for the spin- 
odal. Therefore, we can describe the scaling of events in 
the CA model in the language of cluster scaling for Ising 
models. 

We first discuss how the cluster structure relates to 
thermal critical phenomena in non-mean-field systems 
(R = 1) where hyperscaling is valid. In this case, the 
mean number of clusters in a region of volume £ d , where 
£ is the correlation length, is one. In such systems the 
critical phenomena fluctuation in this volume is isomor- 
phic to the cluster This picture is altered in mean- 
field systems. For mean-field Ising systems (R — > oo) 



there is a line of spinodal critical points in addition to 
the usual critical point. These mean-field thermal singu- 
larities can also be mapped onto percolation transitions 
JlO| | , but the relation between percolation clusters and 
critical fluctuations is qualitatively different. The mean 
number of clusters in a volume £ d is N c = ]{ d e 2 - d / 2 near 
the critical point and N s — R d Ah 3 / 2 ~~ d / 4 near the spin- 
odal Here e = (T-T c )/T c , where T c is the critical 
temperature, and Ah = h — h s , where h s is the value of 
the magnetic field at the spinodal for a fixed temperature 
T < T c . The factor R d appears because all lengths are 
in units of the interaction range. The Ginsburg criterion 
for mean-field critical points is /{R d e 2 ^- dv ) « 1 [§. 
That is, the system is well approximated by mean-field 
theory if the fluctuations are small compared to the or- 
der parameter. Using the mean-field exponents |@] 7 = 1, 
(3=1/2 and v — 1/2, the Ginsburg criterion is equivalent 
to JV C >> 1. We will refer to systems with N c >> 1 but 
finite as near-mean-field. A similar argument is used near 
the spinodal to show that N s >> 1 for near-mean-field. 

Because N c » 1, the meaning of order parameter 
scaling is changed. For systems with hyperscaling [0, 
the density of the single cluster with diameter £ scales as 
e@ , as does the order parameter. In mean- field and near- 
mean-field systems, cannot be the density of a single 
cluster, because that would lead to a magnetization per 
spin greater than one. Instead e 13 is the density of all the 
spins in all the clusters in a volume £ d pT[-p3|. Because 
all of the clusters are identical, the density of each of 
these clusters is ~ e 1 ' 2 /(R d e 2 ~ d ' 2 ) at mean-field crit- 



ical points and p ! s c ~ Ah^ 2 /(R d Ah 3 / 2 -^ 4 ) at spinodals. 
These densities are good approximations in near-mean- 
field systems. We will refer to these clusters as funda- 
mental clusters. These clusters are not the critic al phe - 
nomena fluctuations, but are related to them [prl , ^2|Jl^1 . 

Spinodals mark the boundary between the metastable 
and unstable states. In near-mean-field systems the spin- 
odal is not a sharp singularity but becomes a smeared out 
region [flfjf associated with singularities in complex tem- 
perature and magnetic field space jlffl. As the spinodal 
is approached so is the limit of metastability JlJ]. Hence, 
we would expect that nucleation events, which form an- 
other class of clusters, also play a role in the CA model. 
From the Langevin approach we find that the nu- 
cleation clusters are local regions of growth of the stable 
high stress phase in the metastable low stress phase. An 
earthquake represents the stress release due to the decay 
of the high stress phase into the metastable low stress 
phase. Because the nucleation phenomena of interest oc- 
curs near the spinodal, the classical picture is not valid 
M,19|. Instead, a calculation of the nucleation rate must 
include the effect of the spinodal which involves a vanish- 
ing of the surface tension |20) . With these considerations 
the nucleation rate, which is proportional to the number 
n of clusters per unit volume, is given by p0[| 
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Ah 1 / 2 exp(-AR d Ah 3 / 2 - d / 4 ) 



(1) 



where A is a constant independent of R and Ah. 

The nucleation rate in Eq. ([!]) contains an exponen- 
tial term whose argument is the nucleation barrier. The 
static prefactor, which is independent of the dynamics of 
the model, is l/£ d , where £ = RAh~ x / 4 is the correlation 
length near the spinodal Jl|,|o|,|l) . The Ah 1 / 2 term is 
the kinetic prefactor and is dynamics specific Jl^,^] . For 
the CA model the distance from the spinodal is measured 
by the amount of stress dissipated, i.e., Ah ~ K^/K 
J22J. Finally, the extra factor of R d in the denominator 
reflects the fact that the theory employs a coarse grained 
time scale || , but our simulations use a time scale based 
on plate updates. Because the coarse graining time is 
proportional to the coarse graining volume R d f|§, this 
extra factor is included in the nucleation rate. 

Our assumption is that the CA model behaves like an 
Ising model near the spinodal for mean-field and near- 
mean-field systems. In this limit [p^|-|L5|, fundamental 
clusters and nucleation events, which involve coalescence 
of fundamental clusters fil| , are the only clusters. Be- 
cause it is the decay of the high stress clusters that is the 
"earthquakes" in this model, cluster scaling and earth- 
quake scaling are the same. To understand how this 
point of view provides answers to the questions posed 
in the introduction we discuss the fundamental clusters. 

In mean-field each block fails at the failure threshold 
and fails only once during an earthquake (|||. This be- 
havior is an excellent approximation in near-mean-field 
H . The amount Act of stress transmitted to a site dur- 
ing the failure of a cluster is proportional to the number 
of cluster sites in the interaction volume R d 1 which is 
p i sR d because we are near the spinodal, times the frac- 
tion of stress transmitted from a failed block, which is 
proportional to R~ d . Hence, Act ~ pf during the fail- 
ure of a fundamental cluster, and the mean size of the 
fundamental cluster is s — p^£ d — Ah^ 1 . The num- 
ber of fundamental clusters per unit volume is nf c = 
Rd Ah 3/2-d/4^d = & h 3/2_ Therefore, the density of fun- 
damental clusters with s blocks scales as nf c ~ 1/s 3 / 2 . 

To identify the fundamental clusters we examine the 
stresses on the blocks that make up a cluster of failed 
sites and determine the minimum stress, <7 m i n , of the 
failing sites prior to failure, but after the update of the 
loader plate that triggers the event. We record only those 
events for which CT m ; n is within the window [ct f — Act, ct f ]. 
In Fig. ^ we plot the log of the number of fundamen- 
tal clusters versus log s. For the chosen parameters, 
Aft = 0.01 so that Act ~ 0.01. The slope of 1.53 is con- 
sistent with the theoretical prediction. The mean size 
s = Ah^ 1 = 100 is also consistent with our data. Note 
the lack of data spread. In Figs. |l| and |]the fundamental 
clusters make up only the small s end of the scaling plot, 
but the fundamental clusters comprise ~ 17 x 10 6 out of 
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FIG. 2. Same data as Fig. |lj with the events separated into 
different classes. The triangles represent fundamental clusters 
with scaling exponent of 1.5 ± 0.05. The diamonds are failing 
nucleation clusters with a scaling exponent 1.5±0.05, showing 
the absence of critical slowing down. The open squares are 
arrested nucleation events with a scaling exponent 2.0 ±0.1 
indicating the presence of critical slowing down. There are 
"breakout" events at the largest end of the plot that do not 
exhibit scaling. 



18 x 10 6 events. Hence, a simulation of earthquake faults 
will require huge numbers of events to probe the statistics 
of the interesting and important large event region. 

We now consider the nucleation events and their clus- 
ters. The size of the nucleation events depends on sev- 
eral factors that determine precisely when a given high 
stress nucleation event will stop growing. We will con- 
centrate on two regimes. The first is events near the top 
of the saddle point hill associated with the barrier be- 
tween the stable and metastable state [|I8|J2C|1 . The reason 
we neglect clusters on all scales between the fundamen- 
tal clusters and the saddle point clusters is that Ising 
model studies have found no clusters in this intermedi- 
ate region |l^,^5|. Another aspect of nucleation events 
of this kind in Ising systems is that one must get very 
close to the spinodal to observe critical slowing down be- 
cause the saddle point hill appears to be high but not 
very flat until the system is very close to the spinodal 
|0,E4|. The absence of critical slowing down near the 
saddle point also has been seen in the CA model p3| . 
Hence, there is a class of nucleation events that do not 
quite reach the top of the saddle point hill. As a result, 
random fluctuations lead to the decay of these clusters 
back to the metastable phase. We call these clusters fail- 
ing nucleation events. The probability of these events is 
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characterized by a saddle point calculation without the 
kinetic prefactor. From these considerations and Eq. (|l]), 
the mean number of failing nucleation events per unit 
volume is n fn oc exp(-AR d Ah 3 ^- d ^ 4 )/^ d , the same as 
Eq. (jl]) without the kinetic prefactor. 

To obtain predictions for the scaling regime two results 
are needed. The first is that sf n = A/i 1 / 2 ^, where Ah 1 / 2 
is the density of the nucleating cluster |2(|. Second, be- 
cause the exponential is a rapidly increasing function of 
its argument, as Ah decreases, the probability of a cluster 
increases from almost zero to a relatively large number 
over a very short interval of Ah. The value of Ah where 
this crossover occurs is the limit of metastability [fl6|| . 
For this reason essentially all of the nucleation events 
take place at a fixed value of AR d Ah 3 / 2 ~ d / i = C. As 
for the fundamental clusters the stress transfer to a site 
in a nucleation event is equal to the density of the event. 
For our parameters the density is Ah 1 / 2 = 0.1. Hence we 
identify these nucleating clusters by selecting only those 
events whose cr m i n falls within the window [0.90,0.99]. 
The size of the event is s fn = Ah 1/2 £ d ~ 1000 and the 
number of these events is nf„ oc e /£ . Using the relation 
between R and Ah implied by a fixed value of C, we find 
that nf n scales as l/sj^ 2 . In Fig. || we show a log-log plot 
of rif n versus Sf n . The slope and mean size are consistent 
with our predictions. 

Finally, we consider a second class of nucleation events. 
These are the events that have made it to the top and 
over the saddle point hill and have become arrested dur- 
ing their growth phase, i.e., after growing to some size 
the high stress nucleation region decays back to the low 
stress metastable state. Because the clusters have made 
it to the top of the saddle point hill, this decay of the high 
stress phase is no longer induced by random fluctuations: 
it appears in the Langevin approach due to a decreasing 
loader plate velocity on the coarse grained scale P]22| 
which pulls the system away from the spinodal. We will 
call these clusters arrested nucleation events. 

Because these clusters experience critical slowing 
down, their number per unit volume is given by Eq. (Q). 
A key feature in the growth of nucleation events near 
the spinodal is that their initial growth is a filling in 
p4| , and hence these clusters are compact, that is s an oc 
= R d Ah~ d / A . The density is of order unity so that 
we will identify these events with those clusters whose 
minimum stress of the failing blocks obeys the condition 
Cmin < 0.90. Using the same arguments as in for the 
failing nucleation events, we find that their mean size is 
about 10 4 and the slope of the scaling plot is predicted 
to be 2. The data presented in Fig. || is consistent with 
these predictions. 

In summary, these theoretical considerations and nu- 
merical results strongly suggest several important points. 
(1) Earthquake fault models are statistically dominated 
by small, and in the case of real earthquakes, uninter- 



esting events. (2) The large and small events have dif- 
ferent physical mechanisms. (3) The scaling regime is 
composed of events with two different power law distri- 
butions, which accounts for the data spread at the large 
events end of the scaling plot in Fig. 1. (4) Note that 
there is still a spread in the data at the large events 
end of Fig. ^ and that these events do not scale with a 
slope of 2. Numerical investigation indicates that these 
are "breakout" events that are generated by the spatial 
coalescence of arrested nucleation events [^6| and are be- 
yond the assumptions of our present theoretical treat- 
ment. That is, as the arrested high stress cluster decays, 
it releases stress into the surrounding system. If, due to 
past history, the stress field is unstable, this stress re- 
lease can lead to runaway failure. This type of event was 
considered in Ref. [^5| and is a fourth mechanism that 
must be considered in the generation of earthquakes. In 
contrast, the nucleation events are generated by the co- 
alescence of overlapping fundamental clusters occupying 
the same region of volume £ d . 
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